Frictional dynamics of stiff monolayers: from nucleation dynamics to thermal sliding

Abstract

The inherently nonlinear dynamics of two surfaces as they are driven past each other, a phenomenon known as dry friction, has yet to be fully understood on an atomistic level. New experiments on colloidal monolayers forced over laser-generated substrates now offer the opportunity to investigate friction with single-particle resolution. Here, we use analytical theory and computer simulations to study the effect of thermal fluctuations on the stick-slip mechanism characteristic for the frictional response of a stiff colloidal monolayer on a commensurate substrate. By performing a harmonic expansion of the energy and employing elementary statistical mechanics, we map the motion of the monolayer onto a simple differential equation. Analytical expressions derived from our approach predict a transition from nucleation dynamics, where the monolayer moves in a sequence of activated hops over energy barriers, to “thermal sliding”, in which the effective substrate barrier opposing the motion of the monolayer disappears due to thermal fluctuations, leading to continuous, uninterrupted sliding motion. Furthermore, we find that the average velocity of the monolayer for large driving forces obeys a simple scaling behavior that is consistent with the existence of a static friction. For small forces, however, nucleation provides a mode of motion that leads to a small but non-vanishing mobility of the monolayer. Data obtained from simulations confirm this picture and agree quantitatively with our analytical formulae. The theory developed here holds under general conditions for sufficiently strong inter-particle repulsions and it yields specific predictions that can be tested in experiments.

---

1 Introduction

Although macroscopic laws of friction are centuries old and sufficiently accurate for a multitude of applications,¹ and mesoscopic treatments have yielded great insights into this phenomenon,^2,3 advances in nanotechnology require an understanding of friction on length and time scales in which atomistic details become important.^4,5 Experimentally, several recent approaches, including atomic force microscopy setups,^6–9 quartz-microbalance setups,^10–13 and monolayers of charge-stabilized colloidal particles exposed to light-induced interference patterns,¹⁴ have made it possible to study friction on an atomistic level. In all of these cases, insights are gleaned from the dynamical response of small samples of particles that are driven across an external potential landscape. In particular, the remarkable accuracy and level of control achieved in experiments using two-dimensional crystals in laser fields¹⁴ have shed new light on the dynamics of collective excitations lying at the origin of friction. These experiments have been complemented by simulations and analytical studies of the Frenkel–Kontorova (FK) model.^15–25 Despite its simplicity, this model captures the rich dynamical behavior and, more specifically, the soliton and anti-soliton structures largely determining the magnitude of friction observed in experiments.

Here, our aim is to investigate the effect of thermal fluctuations on the frictional response of stiff monolayers, i.e., monolayers with large elastic moduli, driven over a commensurate substrate. In a previous work,²⁶ we examined the friction induced by a substrate with potential wells arranged in a regular structure identical to that of the unperturbed monolayer driven over it. In particular, we were interested in the role that the inter-particle interaction strength had on the sliding phases that the monolayer adopted as it was driven over such a commensurate substrate. We found that the inter-particle interaction strength regulates how many defects appear in the monolayer as it slides. When the interaction strength is so large that no defects appear in the system, the monolayer adopts a type of stick-slip motion which we call a “hopping wave”. This sliding mechanism consists of long periods of virtual motionlessness which are interrupted by the formation of small nuclei of particles that have escaped from their respective substrate wells. These particles subsequently initiate a cascade of particle hops which eventually encompass the entire system (Fig. 1). The “hopping wave” mechanism is of particular note, since it accounts for a sliding phase in which the monolayer remains structurally intact. It also appears to be a feature of systems of this type because the same mechanism was observed under different conditions in an earlier simulation study conducted by Reguzzoni et al.¹⁷ The purpose of our present work is to provide a general and quantitative understanding of the formation of these hopping waves.

Fig. 1 Snapshot of a stiff monolayer driven by a constant force, F_d, that acts from top to bottom. Each of the particles is colored according to its substrate potential value. Red particles are at the top of their substrate potential barriers, purple particles are at the bottom of the wells, and green particles are somewhere along the walls.

2 Model

We simulated a two-dimensional array of overdamped particles exposed to an external substrate potential while being driven by a constant force. The particles repel each other via Yukawa interactions and receive random kicks due to the solvent surrounding them. The equation of motion for particle i, with position vector r_i, is therefore

The substrate force, F_sub(r_i), acting on particle i is the negative gradient of the externally applied potential, U_sub(r_i) = −(U₀/9){3 + 2[cos(k₁r_i) + cos(k₂r_i) + cos(k₃r_i)]}, where the k-vectors are chosen from the set k_i/‖k‖ ∈ { (0, 1)} with norm . This choice of the k-vectors produces a hexagonal arrangement of potential wells with a lattice constant a and lattice vectors . The depth of the potential minima was set to U₀ = 27k_BT, where k_B is the Boltzmann constant and T is the temperature, and the lattice constant of the substrate was set to a = 6 μm, so that the filling fraction is one, i.e., the number of substrate minima matches the number of particles in the system. As a result of the choice of these parameters, the maximum force that the substrate can exert on a particle in the x direction is F_max = 24πk_BT/a. The colloidal particles we study are charge-stabilized and therefore repel each other via Yukawa interactions. Two particles separated by a distance r from each other have a potential energy of U_Yuk(r) = [capital Gamma, Greek, tilde]e^−κr/r, where the inverse screening length, κ = 6.25 μm⁻¹, determines the range of the Yukawa interaction. The coupling parameter [capital Gamma, Greek, tilde] is related to the effective charge on each colloid. Instead of specifying [capital Gamma, Greek, tilde], we define the interaction strength as Γ = [capital Gamma, Greek, tilde]e^−κa/a, which is the potential energy of two particles separated by one lattice constant. The choice of the parameter values used here resulted in the best agreement between simulation and experiment.^14,26 For particles obeying overdamped Langevin dynamics, at each moment in time a particle experiences an uncorrelated Gaussian random force, Fⁱ_random, with zero mean and variance 2k_BTγ. The friction constant γ is related to the diffusion constant of a single particle by the Einstein relation, γ = D/k_BT. We prepared a rectangular simulation box with periodic boundary conditions compatible with N = 6400 hexagonally arranged substrate minima. At the beginning of each simulation run, we placed a single particle in each well and applied a constant driving force, F_d = (F_d, 0), in the x direction. The algorithm presented in ref. 27 was employed to simulate the motion of the particles with a time discretization of δt = 10⁻⁴. Data were gathered from 100 runs of 2 × 10⁶ time steps for a multitude of values of Γ and F_d while keeping κ, a, and F_max constant. Before performing measurements, the systems were equilibrated for 10⁵ time steps. We used the reduced units γ = k_BT = 1 and all distances were rescaled by the lattice constant a and all forces by F_max.

3 Buildup phase

The stick-slip mechanism mentioned above appears when the inter-particle interaction strength Γ is large enough, in comparison to the substrate potential depth, to preserve the hexagonal structure of the monolayer at all times. The monolayer, if it is to slide, has no other recourse than to form a distortion wave that travels through the system until each particle has moved forwards by one lattice constant. This is accomplished by forming a localized cluster of particles that escape their respective substrate wells. This cluster of hopping particles is the source of a sequence of particle hops that travels through the entire system. After such a “hopping wave” has run its course, the monolayer re-equilibrates and awaits the formation of yet another cluster (see ESI videos V1 and V2†).

In Fig. 2, this process is illustrated for a monolayer driven at two slightly different values of F_d. In the top row of the figure, the average displacement of each particle from its position at t = 0 is plotted as a function of time. One can clearly see that the trajectories alternate between a “buildup phase”, in which the monolayer is in the process of forming a hopping wave, and a “slipping phase”, in which a hopping wave travels through the system until each particle has moved forwards by one lattice constant. Although the applied driving forces differ by only a fraction of a percent, the velocity of the monolayer changes almost by a factor 8. Moreover, the regularity with which the hopping waves appear changes drastically. The plateaus in the graph are colored in blue and correspond to data points belonging to the “buildup phase”.

Fig. 2 Motion of a monolayer with Γ/k_BT = 0.804. We have depicted the absolute displacement of the system, d (first row), the “periodic position of the center of mass”, R/a (second row), the variance in the x direction of the particle positions from the position of the center of mass, σ_x² (third row), and the net force acting on the monolayer, F_net/F_max (bottom row), all as a function of reduced time, γt, for two values of F_d/F_max. Data points colored in blue indicate that the monolayer is in the “buildup phase” whereas data points are colored red when a hopping wave is traveling through the system. Panels in the same column belong to the same trajectories. For animations of the data, see ESI videos V1 and V2.†

We can take advantage of the periodicity of the substrate by considering the “periodic position of the center of mass” R (second row of Fig. 2), which is defined as the average displacement in the x direction of each particle from its nearest substrate potential minimum. During the buildup phase, the value of R is equal to the average displacement of the monolayer, d, modulo the lattice constant a. During the slipping phase, R does not have a physically meaningful value but the sharp peaks are nonetheless clear indicators of the existence of a hopping wave. The value R = 0.25a is of particular significance since it is the point of maximum resistance of the substrate. The plots of R vs. γt not only reveal the repetitive nature of the hopping wave mechanism, but also show that the monolayer on the left gets pinned by the substrate and takes a variable amount of time to evolve a hopping wave, whereas the buildup phase of the monolayer on the right consists of an essentially continuous drift towards the point of maximum resistance of the substrate, and shortly thereafter forms a hopping wave.

One can measure how accurately R describes the positions of the particles in the monolayer by calculating the variance in x direction of the particles' positions, σ_x², from the periodic center of mass (third row of Fig. 2). As can be expected, during the buildup phase, the particles perform small oscillations about R whereas when a hopping wave is traveling through the system, large deviations are observed. We therefore consider the monolayer to be in the buildup phase if the variance is below a cutoff of 8.5 × 10⁻⁴a², which evidently captures the plateaus in the first rows of the columns in Fig. 2.

An examination of the net force acting on the monolayer (bottom row of Fig. 2) confirms that the monolayer driven at F_d/F_max = 0.980 experiences zero net force for significant periods of time, before a random fluctuation creates a hopping wave. The monolayer driven at a rate of F_d/F_max = 0.984, on the other hand, is perpetually in motion. This is surprising since the substrate potential is clearly capable of applying a larger restoring force on each particle than F_d, yet the monolayer glides unhindered in the driving direction. In the following section, we will explain the origin of this “thermal sliding”.

4 Mapping onto the harmonic crystal

The trajectories in Fig. 2 differ primarily in the character and duration of their buildup phases. It also turns out that this part of the trajectories is amenable to analytical treatment. Consider the equation of motion of the periodic center of mass of the monolayer during the buildup phase,where the stochastic force F̄_random(t) is the average random force acting on the particles and is a Gaussian uncorrelated noise with zero mean and variance 2k_BTγ/N. The Yukawa forces cancel due to Newton's third law and so we obtain a stochastic differential equation describing the motion of a Brownian diffuser in a potential defined by the driving force and the average substrate force acting on the monolayer. Since there is no diffusion in the direction perpendicular to the driving force (in the y direction), we merely need to consider the x component of the equation of motion of the system. The only unknown in the equation is the average substrate force acting on the monolayer in the x direction, , which depends on the distribution of the particles' positions r_i and the instantaneous position of the center of mass, R.

In the following, we consider an ensemble of various realizations of configurations with a given R, as they appear in the system at different times along different trajectories. It then follows that the mean velocity of the monolayer during the buildup phase is

where the average 〈…〉_R is performed over multiple buildup phases and F^eff(R) is the effective force that acts on the center mass of the monolayer. To find an analytical expression for F^eff(R), we propose the following two assumptions:

(1) During the buildup phase, the monolayer is in quasi-static thermodynamic equilibrium. It can be expected that this assumption holds because the monolayer creeps along the external potential landscape very slowly during the buildup phase.

(2) The total potential energy of the system can be approximated by a second order Taylor expansion. This approximation is expected to be valid for stiff crystals.

If the first condition is met, then, for a given position R, the probability distribution of the particle positions is the equilibrium distribution with the restriction :

where β = 1/k_BT, δ is the Dirac delta function, and r⃑ is a 2N-dimensional vector of x and y coordinates of each particle. U_tot(r⃑) denotes the total potential energy of the system as a function of all particle positions, and the normalization is given by the 2N-dimensional integral, . The effective substrate force, F^eff(R), defined in eqn (3), is the corresponding ensemble average. The normalization and expectation values of such distributions typically defy analytical evaluation, so a further approximation is required.

The second assumption implies that the position of a particle deviates from the location of the center of mass, R = (R, 0), only by a small displacement, u_i. We can therefore rewrite the position of each particle as r_i = u_i + R + R_i, where R_i is the position vector of the substrate minimum closest to the particle i. The center of mass of the monolayer is able to oscillate in the y direction but, by construction, the deviations in the x direction must cancel, . For small u_i, the total potential of the system can then be approximated by a second order Taylor expansion. Although conceptually simple, the following calculation is rather cumbersome, thus, we have detailed every step of the procedure in the ESI Appendix S1† and mention only the essential results in the following. The work by Baumgartl et al.²⁸ on this approach is recommended for further reading, and an introduction to the harmonic crystal can be found in ref. 29. The second order Taylor expansion of the total potential energy for small u_i is

where the 2N × 2N interaction matrix [Doublestruck D](R) encodes how each particle interacts with the substrate and with all other particles and u⃑^T is the transpose of the 2N-dimensional vector of the small displacements of the particle positions. Both in simulation and in this calculation, nearest neighbor cutoffs were used and, hence, the majority of the entries in [Doublestruck D](R) are zero. The function G(R) combines the linear term in the Taylor expansion of the substrate potential and the driving force F_d, and H(R) is the potential energy of the system if all u_i are zero. Since the interaction matrix is symmetric, it can be brought to a diagonal form. Using the discrete Fourier transformation, the 2N × 2N matrix [Doublestruck D](R) can be resolved into N 2 × 2 matrices.²⁸ The total energy of the system can therefore be rewritten aswhere ū(q) and (q, R) are the discrete Fourier transforms of u and [Doublestruck D](R), respectively, and q is an element of the set of N reciprocal lattice vectors in the Brillouin zone of the substrate lattice. It turns out that the sum over all u_i^x appearing in eqn (5) is , that is, the x component of the Fourier amplitude assigned to q = 0.

Recasting the total potential in this form allows us to evaluate Z(R), the expectation values of the variances of the particles' positions, and the average force that the substrate exerts on the monolayer for a given value of R. The expectation value of the force exerted by the substrate in the x direction is

where ∫du⃑ and ∫dū(q⃑) are 2N-dimensional integrals over all independent degrees of freedom of the system. The only coordinate appearing in the delta function is ū^x(0) and therefore the corresponding integral over this degree of freedom is unity. Since the linear term G(R) is only coupled to ū^x(0), the result of the integration does not depend on G(R). The function H(R) can be pulled out of the integral and cancels with the corresponding term appearing in Z(R). What remains in eqn (7) is the product of 2N − 1 Gaussian integrals that can be evaluated individually for each degree of freedom. Having calculated the variances of u^x, u^y, and their cross correlation, we arrive at the result:where σ_x²(R) and σ_y²(R) are the expectation values of the variances of the particles from their mean position in x and y directions. It turns out that, in the harmonic approximation of the hexagonal lattice, the cross correlation σ_xy vanishes which is compatible with our findings from simulation. These variances are directly related to the diagonal elements of the inverse of the dynamical matrix (q):where μ ∈ {x, y}. The Kronecker delta δ_yμ term arises from the fact that the monolayer is free to oscillate in the y direction and ensures that an additional term is added to the sum when σ_y²(R) is calculated.

We find that the expectation values of the effective substrate force and the variances in the particles' displacements are independent of the applied driving force. Furthermore, the formulae for the variances in eqn (9) can be interpreted as the average value of a function and therefore the effective resistance due to the substrate, F^eff, is an (essentially) intensive quantity. The functional form of F^eff(R) is similar to the original external potential except that it is exponentially reduced in terms of the variances, σ_μ²(R), and thus the temperature T. For a sinusoidal substrate, each of the variances increases monotonically as the monolayer travels along the barrier. We included an instruction file, code for a C program, a python script, and a pair of sample files in the ESI Appendices S3–S7,† with which these formulae can be evaluated.

In Fig. 3, we compare the curves generated by eqn (8) and (9) with data obtained from simulation during the buildup phase. We used 4 monolayers of varying interaction strengths Γ, and applied various driving forces close to, but less than F_max. Data from simulations were compiled by calculating the variances of the particle positions and the net substrate force acting on the center of mass, all as a function of the center of mass R, during the buildup phase. The mean values of these quantities were then plotted along with the theoretical predictions in Fig. 3. A more detailed exposition of how the data were gathered and analyzed can be found in the first section of the ESI Appendix S2.† As predicted by our calculations, there is a range of R/a in which the driving force does not influence the expectation values of the monolayer, which explains the collapse of the data points of the same color within the indicated region. Furthermore, the lines, which are our theoretical predictions, conform very well with the simulation data. We consider the main source of error to be the truncation of the Taylor expansion at the second order because the theoretical curves become more accurate as the interaction strength is increased. The largest error is incurred in the estimation of σ_y²(R) when Γ/k_BT = 0.804, for which the average distance between a data point and the curve is 13% and the corresponding error in the estimation of F^eff(R) is 6.6%. Although the region of space in which our formulae are valid may seem small, the monolayers reside in this region for the vast majority of the time. The harmonic approximation diverges shortly after R = 0.25a because the curvature of the external potential landscape changes sign and therefore ceases to be a pinning potential. This change of curvature explains why, in simulation, hopping waves form almost immediately after the monolayer reaches that point (see right panels of Fig. 2).

Fig. 3 Comparison of theoretical prediction (lines) with data gathered from simulation runs (symbols). Symbols and lines of the same color correspond to the same value of Γ, whereas different symbols of the same color correspond to different driving forces. The dashed black line in the bottom panel of the graph is F^eff(R) for T = 0 or, equivalently, Γ = ∞. Vertical black lines limit the region in which the presented theory is applicable.

5 Dynamical phases

In the previous section, we mapped the buildup phase of the motion of the monolayer onto the motion of a single overdamped Brownian diffuser subject to an effective substrate force, F^eff(R) (eqn (2) and (3)). This effective substrate force is necessarily weaker than the substrate force acting on each particle due to the thermal motion of the particles in the monolayer. The formula for F^eff(R) can be used to find the value of the largest restoring force due to the effective substrate, F^eff_max, which delimits two different dynamical regimes. In the first regime, when F_d > F^eff_max, the monolayer slides unhindered over an effective substrate. Hence we call this motion “thermal sliding”. In the second regime, when the driving force is lower than the effective barrier, F_d < F^eff_max, the monolayer becomes pinned and its center of mass fluctuates about an equilibrium position which can be inferred from Fig. 3. In this latter case, a critical number of particles need to be kicked out of their respective potential wells in order to initiate a hopping wave. The time taken to form such a critical cluster needs to be treated within the framework of the classical nucleation theory.¹⁷

The entire trajectory of a monolayer can therefore be resolved into three phases: the buildup phase, the nucleation phase (where applicable), and the hopping wave phase. We expect that the distributions of both the time taken to complete the buildup phase and the time for a hopping wave to travel through the system are Gaussians. Since the convolution of two Gaussian distributions is also a Gaussian, we define t̃ to be the mean time that the monolayer takes to complete the buildup phase plus the mean time the hopping wave takes to travel through the system, and the quantity ζ² is the sum of the variances of the aforementioned times. The time taken for nucleation to occur, on the other hand, obeys an exponential distribution with a characteristic time τ. The distribution of the total time, t_total, that the monolayer takes to travel forwards by one lattice constant is an exponentially modified Gaussian distribution arising from the convolution of these distributions:

where erfc(x) is the complementary error function and Π is the product of the norms of the distributions. In the top panel of Fig. 4, we have plotted the distribution of the waiting times between two successive hopping waves, P(t_total), for a monolayer driven with a force less than, equal to, and greater than F^eff_max. The data were gathered by finding the time between two successive peaks in the net force acting on the monolayer (see bottom row of Fig. 2). The solid lines are fits of eqn (10) to the data sets, and evidently the tail of the distribution (determined by the value of τ) disappears when F_d becomes larger than F^eff_max. In the bottom of Fig. 4, we plotted the mean nucleation time, τ, gathered from fits of eqn (10), for different monolayers as a function of the driving force. Due to the exponential dependence of τ on the driving force, and the fact that the nucleation barrier is expected to vanish if F_d > F^eff_max, one can observe a change of at least 1 order of magnitude in the nucleation time within a window of 0.4% of F^eff_max. The shape of the curves in the bottom of Fig. 4 indicates that the dynamical transition from nucleation to thermal sliding is continuous.

Fig. 4 Top: distribution of waiting times between hopping waves for a monolayer with Γ/k_BT = 0.804 for driving forces below, equal to, and above the effective barrier height, F^eff_max/F_max = 0.982. The red and black curves have been shifted horizontally by 15γt and 40γt, respectively, for clarity. Bottom: nucleation time, γτ, extracted from the distribution of waiting times for multiple monolayers driven at different rates.

In Fig. 5, we have plotted the mean velocity of a number of monolayers driven at rates both above and below their respective effective force barriers. All of the considerations herein lead to the expectation that there are two different scaling regimes of the mean velocity of the monolayer, γν, with respect to the driving force. In the nucleation regime, also referred to as the creep regime, the value of the mean nucleation time, τ, has the largest influence on the velocity of the monolayer and scales exponentially in the height of the free energy barrier posed by the substrate. The free energy barrier, in turn, is influenced by the driving force and thus the velocity is expected to decay exponentially as F_d goes to zero. In the thermal sliding regime, the mean velocity of the monolayer depends primarily on the time the center of mass takes to diffuse along the potential landscape associated with F^eff(R). If one ignores the random force in eqn (2) and assumes that the variance in the particle positions remains constant, then the first order differential equation for the motion of the monolayer can be solved quite easily. The mean velocity of the monolayer for the thermal sliding regime is then given by γv = [F_d² − F^eff_max²(Γ)]^1/2. The lines in Fig. 5 are plots of this simplification, and conform surprisingly well with the simulation results for F_d > F^eff_max(Γ), especially with respect to the scaling of γν. As a corollary to this consideration, if, for some reason one were unable to measure velocities in the regime close to F_d = F^eff_max, where velocities tend to be very small, then the obtained data would suggest the existence of a static friction located at F_d = F^eff_max(Γ). This apparent static friction obeys Amontons' law, in that the value of F^eff_max(Γ) is independent of the contact area, which in this model is the particle number N, and is, to the first order, proportional to the applied load, which is F_max ∝ U₀.^1,30 Nonetheless, we have shown that the atomic details of such monolayers induce a dramatic change in the scaling of the mean velocity close to F_d = F^eff_max and that it decays exponentially for small driving forces. This finding may have some bearing on the discussion of the molecular origin of static friction found in the literature.^{17,21,30–32}

Fig. 5 Velocity profiles of the simulated monolayers as a function of the reduced driving force F_d/F_max. The lines are plots of the function γv = [F_d² − F^eff_max²(Γ)]^1/2 that terminate at F_d = F^eff_max(Γ) and roughly reproduce the velocity of the monolayers when F_d > F^eff_max(Γ). Inset: mean velocity of a monolayer with Γ/k_BT = 1.038 for different system sizes N.

According to nucleation theory, the mean time of observing a nucleation event is inversely proportional to the product of the nucleation rate and the volume of the system τ = [JV]⁻¹. Therefore, when a collective mechanism is responsible for the sliding motion, one can expect that an increase of the system size will result in an increase of the monolayer mobility, contrary to the traditional experiences with friction. The mean velocity is expected to converge, however, for sufficiently large volumes or high nucleation rates, due to the appearance of multiple, simultaneous, hopping waves. The inset in Fig. 5 is a plot of the mean velocity of a monolayer with Γ/k_BT = 1.038 as a function of F_d for different system sizes N, and, as one can clearly see, the mobility of the monolayer increases with N in the nucleation regime, which is expected to end at F^eff_max/F_max = 0.9855, and indeed, shortly thereafter, the velocities converge. A more detailed examination of the behavior of the system as a function of N can be found in the second section of the ESI Appendix S2.†

We attempted to find a suitable criterion to determine under which conditions the hopping wave mechanism gives way to defect-driven motion, but we have been, so far, unsuccessful. We do expect, however, that the ratio of the substrate potential depth (which favors the formation of defects) and the interaction strength (which penalizes defects) is to the first order the main factor in determining whether the system moves through the formation of hopping waves or through correlated defects. Part of the difficulty in finding a cutoff between these two mechanisms stems from the fact that this is a continuous transition, as illustrated in Video S3,† which shows that monolayers with Γ/k_BT = 0.4 produce both hopping wave and defect induced motion. We summarize our findings in Fig. 6, where we have plotted F^eff_max(Γ) for various values of aF_max/Γ (red dashed line) to delineate the nucleation regime from the thermal sliding regime. The symbols indicate interaction strengths and driving forces that we simulated. Empty boxes represent runs in which the net force acting on the monolayer was zero for significant periods of time, as is necessary for nucleation to occur. Filled circles, on the other hand, denote simulations in which the net force was positive virtually all the time, as expected for thermal sliding.

Fig. 6 Proposed dynamical phase diagram of overdamped monolayers driven by a constant force. The dashed red line is the theoretical demarcation between thermal sliding and nucleation-induced motion, given by F^eff_max(Γ). Empty boxes denote parameter values for which nucleation of hopping waves was found and filled circles correspond to simulation runs in which thermal sliding occurred. The green gradient indicates the region in which the hopping wave mechanism gives way to defect-induced motion.

6 Conclusions and outlook

In this work, we found that stiff overdamped monolayers, driven over commensurate substrates, adopt one of two different dynamic phases and we developed a quantitative theory to describe them. In particular, the velocity of the monolayer can be entirely characterized by the time the center of mass creeps along the effective substrate potential, the time a hopping wave travels through the system, and, where applicable, the mean nucleation time. While the results we presented in this work can be used to determine the first of the aforementioned times, analytical expressions for the hopping wave velocity and the nucleation time need to be developed. With these three quantities, the velocity profiles of the monolayers driven over a commensurate substrate can be reconstructed analytically. The driving force determines which of these times has the largest influence on the velocity of the system and a crossover occurs at F_d = F^eff_max(Γ, F_max, a, β), which has been computed analytically. For large driving forces, the mean velocity of this model scales as if there existed a static friction and, for small forces, we indicated that the mobility has an atypical dependence on the contact area, but a detailed analysis of this finite size effect is still pending.

The theory presented herein remains entirely unchanged if a different radially symmetric interaction potential is employed, provided that the particles in the system always repel each other strongly. We can therefore predict that a density dependent dynamical transition occurs for non-monotonic potentials even for large interaction strengths. The application of the approximation to substrates with different geometries, on the other hand, is much more difficult. It stands to reason, that if the geometry of the substrate is different from that of the monolayer, then the theory is only applicable if, during the buildup phase, each substrate minimum is occupied by exactly one particle. The work conducted by McDermott et al.²² and by Reguzzoni et al.¹⁷ are two examples of systems in which the shapes of the substrate and monolayer are very different, but the theory may remain applicable. Furthermore, although the dynamic phases of the more complex underdamped Langevin dynamics would introduce an additional parameter to the system, it ought to be manageable under the right conditions, particularly in the “onset of sliding” simulations performed in ref. 20.

Acknowledgements

This research was funded by the FWF project numbers P24681-N20, V 305-N27, and P22087-N16. All simulations were performed on the Vienna Scientific Cluster (VSC). The authors gratefully acknowledge Philipp Geiger, Georg Menzl, and Michael Grünwald for useful discussions.

References

1. G. Amontons, Mem. Acad. Roy. Sci., 1699, 206–222.
2. N. Okamoto and M. Nakazawa, Int. J. Numer. Methods Eng., 1979, 14, 337.
3. R. S. Sayles, Tribol. Int., 1996, 29, 639.
4. M. Urbakh and R. Meyer, Nat. Mater., 2010, 9, 8.
5. Y. Mo, K. T. Turner and I. Szlufarska, Nature, 2009, 457, 1116.
6. B. Saha, E. Liu and S. B. Tor, Nanotribological phenomena, principles and mechanisms for MEMS, Springer-Verlag Berlin Heidelberg, 2013.
7. J. A. Ruan and B. Bhushan, J. Tribol., 1994, 116, 378.
8. M. Langer, M. Kisiel, R. Pawlak, F. Pellegrini, G. E. Santoro, R. Buzio, A. Gerbi, G. Balakrishnan, A. Baratoff, E. Tosatti and E. Meyer, Nat. Mater., 2014, 13, 173.
9. R. Pérez, Nat. Mater., 2014, 13, 118.
10. J. Krim, D. H. Solina and R. Chiarello, Phys. Rev. Lett., 1991, 66, 181.
11. T. Coffey and J. Krim, Phys. Rev. Lett., 2005, 95, 076101.
12. L. Bruschi, A. Carlin and G. Mistura, Phys. Rev. Lett., 2002, 88, 046105.
13. J. Jupille, J. J. Ehrhardt, D. Fargues and A. Cassuto, Faraday Discuss. Chem. Soc., 1990, 89, 323.
14. T. Bohlein, J. Mikhael and C. Bechinger, Nat. Mater., 2012, 11, 126.
15. Y. Frenkel and T. Kontorova, Phys. Z. Sowjetunion, 1938, 13, 1.
16. L. Prandtl, Z. Angew. Math. Mech., 1928, 8, 85.
17. M. Reguzzoni, M. Ferrario, S. Zapperi and M. C. Righi, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 1311.
18. O. M. Braun, A. R. Bishop and J. Röder, Phys. Rev. Lett., 1997, 79, 3692.
19. J. Tekić, O. M. Braun and B. Hu, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 71, 026104.
20. O. M. Braun, M. V. Paliy, J. Röder and A. R. Bishop, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 63, 036129.
21. A. Vanossi, N. Manini and E. Tosatti, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 16426.
22. D. McDermott, J. Amelang, C. J. Olson Reichhardt and C. Reichhardt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 062301.
23. C. Reichhardt and C. J. Olson Reichhardt, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2012, 85, 051401.
24. O. M. Braun and Y. S. Kivshar, The Frenkel–Kontorova model: concepts, methods, and applications, Springer-Verlag Berlin Heidelberg, 2004.
25. A. Vanossi and O. M. Braun, J. Phys.: Condens. Matter, 2007, 19, 305017.
26. J. Hasnain, S. Jungblut and C. Dellago, Soft Matter, 2013, 9, 5867.
27. M. P. Allen and D. J. Tildesley, Computer simulation of liquids, Oxford University Press, 1989.
28. H. H. von Grünberg and J. Baumgartl, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 051406.
29. N. W. Ashcroft and N. D. Mermin, Solid state physics, Brooks/Cole, Cengage Learning, 1976.
30. M. H. Müser, L. Wenning and M. O. Robbins, Phys. Rev. Lett., 2001, 86, 1295.
31. M. H. Müser, Proc. Natl. Acad. Sci. U. S. A., 2010, 107, 1257.
32. G. He, M. H. Müser and M. O. Robbins, Science, 1999, 284, 1650.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4nr01790k

---